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ABSTRACT 


We  have  investigated  the  performance  of  a  variety  of  travel  time  computation  methods  for  application  in  highly 
heterogeneous  three-dimensional  (3-D)  velocity  structures,  such  as  those  found  in  Central  Eurasia.  Recent  3-D 
models  of  the  crust  and  upper  mantle  beneath  Central  Eurasia  (e.g.,  Villasenor  et  al.,  2000)  contain  sharp  vertical 
and  horizontal  velocity  contrasts  and  thin  layers  which  complicate  the  accurate  computation  of  travel  times. 

Techniques  based  on  ray-perturbation  theory  or  the  graph  method  (e.g.,  Pulliam  and  Snieder,  1996;  Nolet  and 
Moser,  1993)  do  not  appear  to  work  well  when  the  velocity  anomalies  are  large  relative  to  the  reference  model.  In 
these  cases  the  reference  ray  can  be  very  different  from  the  minimum-time  ray,  and  the  methods  fail  to  converge 
rapidly  and  in  some  cases  do  not  converge  at  all  to  the  global  minimum. 

Ray  shooting  and  finite -difference  algorithms  are  more  computationally  expensive,  but  provide  accurate  results  in 
the  presence  of  large  velocity  contrasts.  In  the  ray-shooting  method,  the  model  is  parameterized  in  terms  of  smooth 
polynomials  in  all  directions.  Ray-shooting  methods  are  normally  implemented  only  in  2-D  (e.g.,  Cerveny  and 
Psenclk,  1988),  and  do  not  consider  propagation  paths  off  the  sagittal  plane.  The  finite  difference  method  computes 
travel  times  for  first-arriving  phases  for  the  entire  model  (e.g.,  Podvin  and  Lecomte,  1991).  Rays  are  obtained 
afterward  by  back-tracing  from  any  point  inside  the  model  (receiver)  to  the  source.  This  method  is  3-D  and  the 
model  is  parameterized  in  terms  of  constant-velocity  cubic  cells.  The  accuracy  of  the  travel  times  is  determined  by 
the  cell  size.  Computer  memory  effectively  limits  the  minimum  cell  size  for  a  given  region  but  this  constraint  can 
be  overcome  by  computing  travel  times  in  2-D. 

Ray  shooting  and  finite  difference  methods  produce  similar  travel  times,  but  differences  can  locally  be  greater  than  1 
s  largely  due  to  the  differences  in  model  parameterization  between  the  two  methods.  We  report  estimates  of 
uncertainties  in  travel  times  and  station  correction  surfaces  between  these  two  methods  and  also  discuss  the  effect  of 
computing  travel  times  in  2-D  rather  than  in  3-D. 

For  IMS  stations  in  Central  Eurasia  we  present  travel  time  correction  surfaces  relative  to  the  1-D  model  akl35 
predicted  by  the  3-D  velocity  model  of  Villasenor  et  al.  (2000).  At  epicentral  distances  less  than  about  10  degrees, 
the  character  of  the  correction  surfaces  is  dominated  by  the  difference  in  crustal  thickness  between  the  3-D  model 
and  akl35.  Beyond  5-10  degrees,  velocity  variations  in  the  uppermost  mantle  become  increasingly  important. 
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OBJECTIVE 


The  objective  of  this  research  is  to  evaluate  methods  for  computing  travel  times  at  regional  distances  for  three 
dimensional  velocity  models,  and  assess  their  usefulness  to  produce  model-predicted  source-specific  station 
corrections  (SSSCs)  for  IMS  stations  in  Central  Eurasia. 
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Introduction 


Locating  low  magnitude  events  with  the  accuracy  required  by  CTBT  monitoring  is  a  challenging  task.  These  events 
are  normally  recorded  only  at  regional  distances,  where  the  effect  of  crustal  and  upper  mantle  sttucture  on  seismic 
ttavel  times  is  important.  One -dimensional  (1-D),  radially-symmetric  models  such  as  iasp91  (Kennett  and  Engdahl, 
1991)  or  akl35  (Kennett  et  al.,  1995)  cannot  predict  travel  times  accurately  for  local  and  regional  distances.  A 
method  to  account  for  the  effects  of  lateral  heterogeneity  at  a  given  station  is  to  consttuct  source-specific  station 
corrections  or  SSSCs.  Empirical  SSSCs  can  be  constructed  in  regions  where  calibration  data  such  as  explosions 
and/or  well  located  earthquakes  are  available  (e.g.,  Bondar  et  al.,  1999).  Unfortunately  this  is  normally  not  the  case, 
and  large  regions  of  the  Earth  have  little  or  no  calibration  data.  For  these  regions,  SSSCs  can  be  obtained  using 
available  regional  and/or  global  velocity  models.  To  obtain  these  model-predicted  SSSCs,  rays  have  to  be  traced 
through  the  model,  from  each  station  to  points  on  a  specified  rectangular  latitude/longitude  grid,  within  a  given 
distance  from  the  station  (normally  20  degrees).  Therefore,  accurate  and  efficient  methods  to  calculate  travel  times 
in  3-D  velocity  models  are  required. 

Computation  of  travel  times  in  heterogeneous  media  is  a  complicated  problem.  This  complexity  is  illustrated  by  the 
large  number  of  methodologies  developed,  and  by  the  fact  that  each  one  has  its  advantages  and  shortcomings  (for  a 
review  see  for  example  Cerveny,  1987).  We  have  restricted  this  study  to  three  families  of  methods:  ray  bending,  ray 
shooting,  and  finite  difference  methods.  Bending  methods  are  generally  fast  and  work  well  for  smoothly  varying 
velocity  models,  but  may  converge  toward  local  travel  time  minima.  Shooting  methods  work  better  for  more 
heterogeneous  models,  but  at  the  expense  of  computation  time,  and  they  cannot  ttace  rays  in  shadow  zones.  Finite 
difference  methods  can  provide  very  accurate  results  even  in  presence  of  extreme  velocity  variations,  but  require 
small  cell  sizes.  Finite  difference  methods  normally  produce  ttavel  times  for  first  arriving  phases,  and  not  for  later 
phases. 

The  two  bending  codes  tested  in  this  study  are  described  by  Bijwaard  and  Spakman  (1999).  The  first  code  is  an 
implementation  of  the  graph  method  by  Nolet  and  Moser  (1993)  based  on  the  work  of  Moser  (1991).  The  second 
code  is  based  on  ray  perturbation  theory,  and  was  developed  by  Pulliam  and  Snieder  (1996)  based  on  the  original 
implementation  of  Snieder  and  Sambridge  (1993).  Both  codes  perform  3-D  ray  tracing  and  work  well  for  models 
with  small  velocity  perturbations  (e.g.  Bijwaard  et  al.,  1998).  However,  they  do  not  work  with  other  models  (e.g., 
Villasenor  et  al.,  2000)  that  contain  sharp  vertical  discontinuities,  thin  layers,  and  large  lateral  velocity  variations. 
These  two  codes  have  not  been  used  in  the  results  presented  in  this  study. 

We  have  also  tested  the  ray  shooting  code  of  Cerveny  and  Psencik  (1988).  This  is  a  2-D  method,  for  laterally- 
varying  layered  models.  The  geometry  of  the  layers  is  described  in  terms  of  splines.  Computations  are  made  in  a 
cartesian  (flat  Earth)  reference  system,  and  applying  an  earth-flattening  transformation. 

Finally  we  have  tested  the  finite  difference  technique  of  Podvin  and  Lecomte  (1991).  This  method  can  be  used  both 
in  3-D  and  2-D.  The  model  is  parameterized  in  terms  of  cubic  constant-velocity  cells.  Therefore,  models  described 
using  a  spherical  parameterization  have  to  be  mapped  into  a  cartesian  reference  system,  resulting  in  some 
inefficiency  in  the  amount  of  memory  required  to  store  the  model.  Travel  time  wavefronts  and  not  rays  are 
computed  in  this  method.  Rays  can  be  obtained  a  posteriori  by  backtracking  towards  the  source  following  the 
normal  direction  to  the  wavefronts. 

Accuracy  tests 

We  have  conducted  tests  to  evaluate  the  performance  and  accuracy  of  the  ray  shooting  and  finite  difference 
methods.  First  we  test  these  methods  against  theoretical  ttavel  times  for  the  one-dimensional,  radially-symmetric 
model  akl35  (Kennett  et  al.,  1995).  Second  we  compare  the  results  of  both  methods  in  2-D,  using  vertical  cross- 
sections  of  the  S  model  of  Villasenor  et  al.  (2000),  and  finally  we  compare  the  performance  of  2-D  versus  full  3-D 
ttavel  times. 

To  evaluate  the  accuracy  of  the  ray  shooting  and  finite  difference  methods,  we  compare  their  results  with  theoretical 
ttavel  times  (Buland  and  Chapman,  1983)  for  akl35.  The  results  of  these  comparisons  are  shown  in  Figure  1.  In 


this  case  the  ray  shooting  method  is  very  precise,  better  that  0.1  s.  The  accuracy  of  the  finite  difference  method 
depends  on  the  cell  size.  For  a  cell  size  of  1  x  1  km  the  precision  is  also  better  than  0.1  s.  If  the  cell  size  is 
increased  to  2  x  2  km  the  precision  is  still  better  that  0.25  second,  but  some  artifacts  caused  by  the  finite  difference 
model  parameterization  (2x2  km  constant  velocity  cells)  can  be  seen.  When  the  first  arrival  changes  from  Pg  to  Pn 
(at  approximately  1.5  degrees)  there  is  a  baseline  shift  of  0.1  s  (Figure  1).  This  shift  is  caused  because  Moho  depth 
in  the  finite  difference  model  parameterization  of  2  x  2  km  cells  is  36  km,  instead  of  35  km  in  akl35. 

We  have  compared  the  2-D  ray  shooting  method  and  the  finite  difference  method  using  two  great-circle  vertical 
cross  sections  through  the  S  model  of  Villasenor  et  al.  (2000)  for  Central  Eurasia.  This  model  is  transversely 
isotropic  in  the  uppermost  mantle,  but  for  this  test  we  have  used  the  vsv  component  of  the  model,  which  is  better 
constrained  than  vSh-  The  differences  between  travel  times  using  both  methods  for  the  two  great  circle  paths 
considered  are  shown  in  Figure  2.  The  differences  are  on  the  order  of  2  seconds,  much  larger  than  the  tests 
performed  with  the  1-D  model  akl35.  We  find  that  these  differences  are  caused  by  the  way  the  3-D  model  is 
translated  into  the  internal  parameterization  of  each  travel  time  code.  In  some  other  cases,  not  shown  here,  the 
differences  are  much  larger,  and  caused  by  the  failure  of  the  ray-shooting  method  to  find  the  minimum  travel  time 
ray  path.  In  most  cases,  modifying  the  parameters  that  control  the  accuracy  of  the  ray  shooting  solution  eliminates 
these  problems,  but  at  the  expense  of  longer  computation  time.  On  the  other  hand,  the  finite  difference  method  is 
very  fast  and  efficient,  and  always  provides  the  first  arrival  time  even  in  presence  of  extreme  velocity  contrasts 
(Podvin  and  Lecomte,  1991). 

Finally  we  have  investigated  the  importance  of  full  3-D  travel  time  computation  versus  computations  in  2-D.  For 
this  test  we  only  use  the  finite -difference  method  of  Podvin  and  Lecomte  (1991)  because  the  ray  shooting  code  is 
strictly  2-D.  To  evaluate  the  differences  between  3-D  and  2-D  travel  times  we  have  considered  two  profiles  through 
the  S  model  of  Villasenor  et  al.  (2000).  Figure  3  shows  the  differences  between  the  computed  3-D  and  2-D  travel 
times.  These  differences  are  small  (less  than  2  seconds),  and  similar  in  magnitude  to  the  differences  between  both  2- 
D  codes  shown  in  Figure  2.  Full  3-D  ray-tracing  is  important  for  non-linear  tomography  studies,  allowing  to  obtain 
improved,  more  focused  velocity  models  (Bijwaard  and  Spakman,  2000;  Widiyantoro  et  al.,  2000),  but  for  the  model 
considered  here  and  for  regional  distances  it  is  not  necessary  for  computing  travel  times  and  SSSCs  with  sufficient 
precision.  This  is  probably  because  the  model  we  have  used  is  more  heterogeneous  vertically  (with  thin  layers  and 
sharp  velocity  contrasts)  than  laterally,  and  therefore  off-great  circle  propagation  is  not  very  important. 


Effects  of  the  reference  model  on  the  computation  of  SSSCs 

Global  velocity  models  obtained  by  inversion  of  large  datasets  of  arrival  times  of  seismic  phases  (e.g.,  Bijwaard  et 
al.,  1988),  are  normally  described  as  perturbations  to  1-D  models  such  as  iasp91  or  akl35.  These  3-D  models 
provide  high  resolution  in  regions  with  high  density  of  earthquakes  and  stations  but  contain  large  gaps  in  regions 
where  seismicity  is  low  and/or  stations  do  not  exist.  On  the  other  hand,  models  that  use  surface  wave  and  normal 
mode  data  have  lower  resolution  than  models  based  on  body  waves,  but  their  coverage  is  more  homogeneous. 
These  models  are  normally  described  as  perturbations  relative  to  PREM  (Dziewonski  and  Anderson,  1981) 

There  are  some  significant  differences  between  PREM  and  akl35  in  the  crust  and  upper  mantle.  PREM  is  designed 
as  an  average  Earth  model,  and  contains  a  3  km  thick  water  layer,  with  Moho  located  at  24.4  km  depth.  The  crust  in 
akl35  is  of  continental  type  (because  most  of  seismic  stations  are  located  on  continents),  with  a  Moho  depth  of  35 
km.  In  the  uppermost  mantle,  between  Moho  and  220  km,  PREM  is  transversely  isotropic,  and  exhibits  a  well- 
developed  low  velocity  zone,  while  akl35  is  isotropic  and  lacks  a  low  velocity  zone,  although  the  velocity  gradient 
in  the  uppermost  mantle  is  very  small  for  both  P  and  S.  In  the  upper  mantle  PREM  has  a  velocity  discontinuity  at 
220  km  depth  which  is  much  sharper  for  P  than  for  S.  This  discontinuity  is  absent  in  akl35.  Finally  the  transition 
zone  is  20  km  thicker  for  PREM  (400  -  670  km)  than  for  akl35  (410  -  660  km). 

These  differences  between  PREM  and  akl35  obviously  translate  into  different  travel  times  at  regional  distances,  as 
shown  in  Figure  4.  We  have  calculated  the  difference  in  travel  times  between  PREM  and  akl35  for  a  source  located 
at  40  km  depth,  and  for  receivers  also  at  40  km  depth.  This  has  been  done  to  avoid  additional  complications  by 
differences  in  crustal  structure.  Therefore,  the  curves  shown  in  Figure  4  represent  the  effect  of  differences  in  upper 
mantle  and  transition  zone  structure  between  the  two  models.  Immediately  below  Moho  vPH  for  PREM  is  larger  than 
vP  for  akl35.  As  a  result  PREM  travel  times  for  vPH  are  faster  (smaller)  than  akl35,  and  this  can  be  observed  as 


increasingly  negative  values  in  Figure  4  (top  panel,  solid  line).  The  difference  in  travel  time  reaches  a  minimum  of 
approximately  -3  s  at  distances  of  14  degrees,  and  then  decreases  to  0  at  25  degrees.  The  situation  is  reversed  for 
PREM  vPV.  Immediately  below  Moho  vPV  for  PREM  is  slightly  smaller  vP  for  akl35  and  this  results  in  positive 
values,  indicating  that  travel  times  for  PREM  vPV  are  slower  (larger)  than  akl35  (top  panel,  dashed  line).  The  same 
pattern  can  be  observed  for  S  wave  travel  times  (Figure  4,  bottom).  In  this  case  the  differences  in  travel  times  are 
larger  (-8  s  for  PREM  vsh)  and  the  maximum  values  of  the  differences  occur  at  different  distances  (e.g.,  at  19 
degrees  for  PREM  vSh)- 

SSSCs  are  normally  referred  to  iasp91  or  akl35  (for  the  purpose  of  regional  travel  times  both  models  are  identical) 
but  many  global  and  regional  models  are  expressed  as  perturbations  relative  to  PREM.  Because  of  the  differences  in 
travel  times  between  PREM  and  akl35  we  must  expect  artifacts  in  the  SSSCs  that  are  caused  by  the  choice  of  the 
reference  model,  in  addition  to  features  due  to  lateral  heterogeneity.  Figure  5a  shows  the  SSSC  relative  to  akl35  for 
auxiliary  IMS  station  AAK  (Ala-Archa,  Kyrgyzstan)  computed  from  the  S  model  of  Villasenor  et  al.  (2000).  At 
epicentral  distances  less  than  about  10  degrees,  the  character  of  the  correction  surface  is  dominated  by  the  difference 
in  crustal  thickness  between  the  3-D  model  and  akl35.  The  crust  in  Central  Eurasia  is  thicker  than  35  km,  and  this 
results  in  positive  station  corrections  (slower  travel  times  than  akl35).  Beyond  5-10  degrees,  velocity  variations  in 
the  uppermost  mantle  become  increasingly  important.  High  velocities  in  the  Indian  craton  produce  large  negative 
corrections  (faster  travel  times  than  akl35).  This  effect  can  also  be  observed  in  the  East  European  platform,  north  of 
the  Caspian  Sea.  On  top  of  these  features  caused  by  lateral  heterogeneity  there  is  a  ring  of  negative  values,  centered 
on  the  station  and  with  a  radius  of  approximately  20  degrees.  As  shown  in  Figure  4,  this  ring  can  be  explained  by 
the  differences  between  our  upper  mantle  model  (which  uses  PREM  as  reference  model)  and  akl35. 

To  avoid  these  artifacts,  the  SSSCs  can  be  calculated  relative  to  the  model  beneath  the  station,  instead  of  relative  to 
iasp91  or  akl35.  This  guarantees  that  the  velocity  discontinuities  of  the  3-D  model  and  reference  model  are  located 
exactly  at  the  same  depths.  Figure  5b  shows  the  results  of  applying  this  procedure  to  AAK.  The  ring  with  negative 
values  has  disappeared,  the  shape  of  the  correction  surface  is  smoother  and  better  correlated  with  the  features  of  the 
3-D  model.  Figure  5c  shows  the  difference  between  the  two  correction  surfaces  obtained  for  AAK  (Figures  5a  and 
5b).  This  clearly  illustrates  that  the  ring  is  an  artifact  caused  by  the  choice  of  the  reference  model. 


CONCLUSIONS  AND  RECOMMENDATIONS 


We  have  tested  the  performance  of  travel  time  codes  based  on  ray  bending,  ray  shooting  and  finite  differences  for 
computing  model-predicted  SSSCs.  For  three-dimensional  Earth  models  with  large  lateral  and  vertical  velocity 
variations  ray  shooting  and  finite  differences  provide  the  best  results. 

For  laterally  heterogeneous  models,  the  differences  in  computed  travel  times  between  different  codes  can  be  rather 
large  (up  to  2  seconds  for  regional  distances).  In  some  cases  this  can  be  caused  by  poor  performance  of  one  the 
methods  (for  example  not  being  able  to  find  the  minimum  travel  time  path,  and  converging  to  a  local  minimum). 
However,  in  the  examples  presented  here,  the  differences  are  caused  by  differences  in  the  internal  model 
parameterization  of  each  travel  time  code.  These  differences  can  be  even  larger  if  the  influence  of  model 
uncertainties  on  the  computation  of  travel  times  is  considered. 

Full  3-D  ray  tracing  is  important  for  non-linear  tomography  studies,  allowing  to  obtain  improved,  more  focused 
velocity  models.  However,  3-D  travel  times  are  not  necessary  for  computing  regional  SSSCs.  Tests  indicate  that, 
for  models  like  the  ones  used  in  this  study,  travel  time  differences  between  2-D  and  3-D  ray  tracing  are  similar  in 
magnitude  to  differences  between  travel  time  codes.  In  the  future,  when  higher-resolution,  more  heterogeneous 
models  become  available,  3-D  ray  tracing  may  be  required. 

The  choice  of  the  reference  model  is  extremely  important  for  determining  the  features  of  model-predicted  SSSCs. 
To  ensure  smoothly  varying  SSSCs,  the  3-D  model  should  have  the  upper  mantle  and  transition  zone  discontinuities 
at  the  same  depth  as  the  reference  model.  Otherwise,  the  computed  SSSCs  may  display  large  artifacts,  that  manifest 
as  circular  features  (positive  or  negative  in  sign)  centered  on  the  station.  These  artifacts  can  be  eliminated  by  using 
the  model  beneath  the  station  as  the  reference  model  for  the  computation  of  the  SSSCs. 


Finally,  methods  need  to  be  developed  to  evaluate  the  uncertainties  of  travel  times  and  SSSCs  based  upon  the 
uncertainties  in  the  travel  time  methods  and  in  the  3-D  velocity  models. 
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Comparison  of  travel  time  methods  for  1-D  model  ak135 
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Figure  1.  Differences  between  computed  travel  times  using  different  methods  and  theoretical  arrival  times  for 
the  one  dimensional  model  akl35.  Thick  black  line  for  the  finite  difference  method  with  1  x  1  km  cells,  thick 
gray  line  for  the  same  method  with  2x2  km  cells,  and  dashed  line  for  the  ray-shooting  method. 


2-D  ray  shooting  vs.  2-D  finite  differences 
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Figure  2.  Differences  between  travel  times  computed  using  the  2-D  ray  shooting  and  the  2-D  finite  difference 
methods  for  two  great  circle  cross  sections  through  the  vSV  model  of  Villaseasor  et  al.  (2000).  The  starting  point 
of  both  great  circle  cross  sections  is  IMS  station  AAK  (Ala-Archa,  Kyrgyzstan).  The  curves  represent  the  ray¬ 
shooting  travel  time  minus  the  finite  difference  travel  time.  Map  inset  shows  the  location  of  the  two  cross 
sections. 


2-D  finite  differences  vs.  3-D  finite  differences 
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Figure  3.  Differences  between  computed  travel  times  along  two  great  circle  cross  sections  through  the  vSV 
model  of  Villaseaeor  et  al.  (2000)  using  the  finite  difference  method  in  2-D  and  3-D.  The  starting  point  of  both 
cross  sections  is  IMS  station  AAK  (Ala-Archa,  Kyrgyzstan).  Each  curve  represents  the  2-D  travel  time  minus  the 
3-D  travel  time.  Map  inset  shows  the  location  of  the  two  cross  sections. 
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Figure  4.  Travel  time  differences  between  1-D  models  PREM  and  akl35  for  regional  distances.  Top  panel 
shows  differences  in  P-wave  traveltimes,  and  bottom  panel  shows  differences  in  S-wave  travel  times. 


Figure  5.  SSSCs  lor  [MS  station  A AK  ( Ala-Archa.  Kyrgyzstan I  For  Ihc  vSV  model  ol'Villasenor  d  al.  12000).  a) 
correction  surface  relative  lo  akIJS.  bi  relative  lo  Ihc  model  beneath  the  station,  and  cl  dilference  between  the 
two  SSSCs. 


